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INTRODUCTION 

The present work on calculation of free shear flows was carried out to investi- 
gate the usefulness of several concepts which were previously successfully applied to 
wall flows. The method belongs to the class of differential approaches. The turbu- 
lence is taken into account by the introduction of one additional partial differential 
equation, the transport equation for the turbulent shear stress. The structure of tur- 
bulence is modeled after Bradshaw et al. (ref. 1). This model has been used success- 
fully in boundary layers and its applicability to other flows is demonstrated in this con- 
tribution. An earlier attempt to use this approach for calculation of free flows was 
made by Laster (ref. 2). The work reported here differs substantially from that of 
Laster in several ways. The most important difference is that the region around the 
center line is treated by invoking the interaction hypothesis (ref. 3) (concerning the 
structure of turbulence in the regions separated by the velocity extrema). The com- 
pressibility effects on shear layer spreading at low and moderate Mach numbers were 
investigated. In the absence of detailed experiments in free flows, the evidence from 
boundary layers that at low Mach numbers the structure of turbulence is unaffected by 
the compressibility was relied on. The present model was tested over a range of self- 
preserving and developing flows including pressure gradients using identical empirical 
input. The dependence of the structure of turbulence on the spreading rate of the 
shear layer d6/dx was established. 


SYMBOLS 

aj,G,L defined by equations (2) 

Cp specific heat at constant pressure 

L Q width to the half velocity point on the profile 

M Mach number 

p mean pressure 

jjt " " — 
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recovery factor 


q2 = u: u : 
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T temperature 

XJ,V mean velocity components 

Ul,U 2 external velocities at edges of a mixing layer 


AU maximum velocity difference across shear layer 

u',v\w' fluctuating velocity components 

~ p'v* 

V = V + 

coordinate axes 
distances along axes 
angle of characteristic 
ratio of specific heats 

V 

shear layer thickness defined as distance between points where r = 0.05r m 

d6/dx 

^ (d6/dx) still _ air j et 

e dissipation rate 

0 momentum thickness 

p mean density 

o spreading parameter for free shear layers 
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r + and t~ shear stress profiles of "simple" layers 

Superscripts: 

a exponent 

' fluctuating quantities 

Subscripts: 

<L center- line value 

i index 

J initial jet value 

m,max maximum value 

1/2 half velocity point 


ANALYSIS 


Equations and the Model of Turbulence 

The governing equations considered are the continuity, momentum, and turbulent 
kinetic energy equations: 
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The turbulent kinetic energy equation contains three additional correlations not 
appearing in the other two equations. To close this system, assumptions about the struc- 
ture of turbulence would have to be made, where by structure a given relation between the 
local values of two turbulent quantities is understood. As in Bradshaw et al. (ref. 1), 
three relations are used to define the structure: 
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(2a) 


r = a x q 2 



(2b) 
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The functions a^, L, and G are specified by algebraic expressions and the local 
length and velocity scales of turbulence are assumed to be proportional, respectively, to 
6 and |r m |l/ 2 . 

By utilizing relations (2), the conservation equation for q2 may be converted into 
an empirical transport equation for t 
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The choice of convective type diffusion (eq. (2c)) over gradient type in free flows is 
supported by mixing layer experiments. (See refs. 4 and 5.) The positions of zero diffu- 
sion and maximum kinetic energy do not coincide, a fact for which gradient diffusion can- 
not account for. The two points are separated by a distance of the order of 5 percent of 
the shear layer thickness. Also, the free shear flows exhibit strong large-scale motions 
which made the convective diffusion important. 

There is much experimental evidence to support the relation (2a). An examination 
of a wide range of experimental data showed that the value of a^ varies within a small 
range depending on the flow considered. The only difficulty in this formulation occurs in 
the vicinity of a velocity extremum, as discussed in the next paragraph. Relation (2b) is 
a logical extension, based on equation (2a), of 


e 



which is a commonly accepted model. 


Flows With Velocity Extrema 

The formulation presented was applied to free mixing layers with the simple choice 
of a^(y) = Constant as in boundary layers . In flows with velocity extrema, the shear 
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stress changes sign in the vicinity of the velocity extremum . The turbulent kinetic 
energy, by definition, does not change sign and this fact precludes the use of constant 
aj. Also, the shear stress equation is singular at the point where a^ = 0 (its solution 
being regular) and for this reason presents numerical difficulties. (See ref. 2.) 

To avoid this problem, the suggestion of Bradshaw (ref. 3) is used and the flow with 
velocity extremum is regarded as two adjoining "simple" shear layers which interact only 
through the mean velocity profile. Each layer has its own shear profile and the algebraic 
sum of these profiles in the region of overlap gives the shear profile of the complete flow. 
One reason for looking at the flow from this viewpoint is that if the structure of turbulence 
in each layer is unaffected by the presence of the adjoining layer (does not actively partic- 
ipate in the interaction) , then a simple tool for calculating more complex flows such as 
jets, wakes, and wall jets is obtained. There are two such simple or basic shear layer 
flows — the boundary layer and the mixing layer - and it is proposed to regard all other 
thin shear layer flows as combinations of these two. The empirical functions in each 
layer of a complex flow were to be the same, or nearly the same, as in the corresponding 
simple shear layer and, as a result, the task of determining them would be simplified. 

(The actual difference between jets and mixing layers can be seen in figures 1 and 2.) 
Another reason is that this point of view allows a simple explanation of the regions of 
"negative production" of turbulent kinetic energy which occur in asymmetric flows. The 
technique can be utilized for calculating these flows which otherwise require the a priori 
knowledge of the point of vanishing shear when the original formulation is used. The 
same shortcoming affects ail models which involve the eddy-viscosity concept in one form 
or another. 

The idea was applied to the duct flows (ref. 3) and very good results were obtained 
with the empirical functions which were developed for boundary layers. In this work the 
same approach is applied to jets and wakes. 

Empirical Functions of the Structure 

The empirical functions of the structure of turbulence in mixing layers (with 
U 2 /U 1 = 0 ) were derived from the experimental data (refs. 4 and 6) and then refined by 
comparison of the velocity and shear profiles with the experiments. (See fig. 1.) The 
simple choice of constant a^ = 0.15 and L/6 = 0.09 was found to be adequate for good 
results. The shape of the diffusion function G was obtained by integration of the diffu- 
sion in the turbulent kinetic energy balance. The proper values of the empirical func- 
tions for U 2 /U 1 * 0 were deduced from calculations by comparison with the results of 
Spencer (ref. 7). The empirical functions were found to be dependent on the velocity ratio 
of the mixing layer. This dependence may be correlated with dS/dx, the spreading rate 
of the shear layer. 
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The same model was applied to jets and wakes as outlined. The two adjoining layers 
are calculated as two separate layers. They share joint U and V profiles but they 
have separate shear profiles. The only difference from the mixing- layer program appears 
in the momentum equation where the shear profiles are superposed to give the "true" 
shear profile in the region of overlap: 
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Beyond the velocity maximum, the layer experiences a negative production in the shear 
equation which limits the region of overlap. 

The program is written, at the present time, for symmetric jets and wakes because 
of lack of experimental data. However, the concept is not restricted in any way to sym- 
metric flows and its importance lies in its ability to treat and explain asymmetric flows. 

The calculations show that in free flows the structure of the "simple" layer is 
affected by the interaction. The interaction tends to modify more the magnitude than 
the shape of the empirical functions. Thus, the and L/6 were retained constant 
and the shape of the diffusion function G was slightly altered. (See fig. 2.) The com- 
parison with experiments indicates again the dependence of the structure on the spreading 
rate d6/dx. The limiting case for d6/dx = 0 is shown in broken lines in figure 2. 

The results show the usefulness of the interaction concept for calculation of free 
shear flows. The required modifications of the empirical functions are not large and 
indicate that the flows with velocity extrema may be regarded as weakly interacting 
adjoining shear layers . There appears to be more interaction of turbulent structure in 
the free flows than in the ducts and the explanation may be sought in the behavior of the 
shear stress in the vicinity of the velocity extremum. (See fig. 3.) The shear stress 
profiles of free flows overlap significantly more and the value of the shear on the center 
line is typically 0.55T max , against 0.1r max in the duct, and thus causes a stronger 
interaction. 

All the results presented are calculated with the same input, the structure being a 
function of y/ 6 and d6/dx alone. 


Compressible Flow 

The governing equations for the compressible flow are much more complicated than 
their incompressible counterpart. H the restriction is made to include only low and mod- 
erate Mach number flows, many of the new correlations appearing in the compressible 
equations may be neglected on the basis of order-of-magnitude arguments. (See ref. 8.) 
This neglect simplifies the equations which may then be written in boundary- layer form 

as follows: 
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The inclusion of compressibility presents two additional tasks . The structure of turbu- 
lence in compressible flow and the density variation across the layer have to be deter- 
mined. 

For the values of the empirical functions in compressible flows, the suggestion 
(ref. 9) is relied on that the turbulence is convected passively by the mean^flow as long 

as the local Mach number of the root-mean-square fluctuations ^pu’^/ypj is much 
less than unity. The inference is that the structure of turbulence is Mach number inde- 
pendent and the ’’incompressible" values of the empirical functions may be used. This 
assumption is based on the analysis of experiments in boundary layers, where this condi- 
tion is satisfied up to moderate Mach numbers (M < 5) . 

In the Mach number range to which the present approach is limited, a good approxi- 
mation for the density profile can be obtained from the Crocco formula used in boundary 
layers: 
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CpT + 0.5rU^ = Constant 

which together with the equation of state yields p = p(U). It eliminates the need for an 
additional equation and may be incorporated into the incompressible version of the pro- 
gram with only small modifications. 

Compressibility effects due to temperature or density differences were also incor- 
porated by assuming, respectively, similarity between temperature and velocity profiles 
and between mass fractions and velocity profiles. That again yields a relation p = p(U) 
which is useful for small compressibility effects. Large effects will require the solution 
of a separate equation. 


METHOD OF SOLUTION 

With the diffusion term in the kinetic energy equation modeled to be of the convec- 
tive type, the set of equations becomes hyperbolic. There are two choices of solving this 
system, either by a procedure suitable for parabolic equations or to use the method of 
characteristics. The first approach is advantageous if it is intended to introduce addi- 
tional equations, for example, for compressible flows or for a more involved model of 
turbulence. In the present work, it was decided to use the mathematically simpler method 
of characteristics used already in the boundary- layer calculations. (See ref. 1.) The 
model of turbulence and the empirical functions developed in the course of this work 
would be the same if an alternative method of solution of these equations is used. 

The accuracy of the calculations is governed by the number of points on the profile 
and by the number of iterations used to improve the interpolation along the characteris- 
tics. However, even if several iterations are used, the momentum and mass balance are 
not preserved because of numerical inaccuracies and lead to a momentum thickness gain 
of the order of 0.1 percent per station. To avoid accumulation of error, the grid size is 
being readjusted by that very small amount after every step so that momentum is pre- 
served. The mass balance is then very well preserved - to within 1 percent on a typical 
run. The largest inaccuracies occur for flows issuing into a small external stream. One 
of the characteristic angles near the edge tends to a= arc tan (V/U) and precludes cal- 
culation of mixing with still air unless some special numerical treatment for this bound- 
ary is introduced. The problem at the still-air edge affects other methods of solution as 
well. Calculation of mixing with small external stream has some practical limitations. 
The step size in the x- direction is inversely proportional to tan and, therefore, 

for economical calculation, it is necessary to maintain at least Uj = 0.05U max which 
gives tan a max ~ 1. The region of very small external stream also suffers from larger 
than average errors in mass balance. 
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COMPARISON WITH EXPERIMENTS 


The model developed for the mixing layers can be used for all velocity ratios. Cal- 
culations for three ratios are compared with experiments below: 

(a) Mixing into still air. (See Liepmann and Laufer (ref. 6) and Bradshaw and 
Ferriss (ref. 4) and also figure 4.) Note that the computed profile has a small external 
stream rather than still air on its edge; thus, the spreading rate is reduced somewhat. 

(b) Mixing with a parallel moving stream. (See Spencer (ref. 7) and figures 5 and 6.) 

The jet program is capable of handling both jets and wakes in the presence of pres- 
sure gradients by using the same empirical input. Some predictions compared with exper- 
iments are 

(a) Jet mixing with still air. (See Bradbury (ref. 10) and figure 7.) Both the exper- 
iment and the calculations had actually a small external stream which does not affect the 
nondimensional profiles . 

(b) Small deficit wake. (See Townsend (ref. 11).) The region in which Townsend 
made his measurements is far from complete self-preservation as documented by the dif- 
ference between the measured shear profile and the profile required for self-preservation. 
(See fig. 8.) The calculations compare well in the region investigated by Townsend and 
show self-preservation far downstream. (See fig. 9.) The small excess jet tends toward 
the same results as the small deficit wake when W — 0. 

When the external flow varies as Uj x a , there are ranges of the exponent a for 
which the momentum equation allows self-preservation of jets or wakes. The following 
two flows fall in that category: 

(c) Wakes in a pressure gradient - investigated by Gartshore (ref. 12) who found 
them approximately self-preserving. The calculations show that, although the flow has 
the tendency to conform to the self-preservation, it drifts steadily away from it (fig. 10). 

(d) Self -preserving jet in a pressure gradient. Value of the exponent a = 
was used for the calculations. The results obtained are plausible, although there are 
no experiments to support them. (See fig. 11.) 
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APPENDIX 


COMMENTS ON TEST CASES 

This appendix presents comments on some of the test cases presented in figures 12 

to 19. 


Test Case 5 (Hill and Page) 

The free shear layer was measured in the early stages of its development from a 
boundary layer separated over a cavity. At the separation the boundary layer was turbu- 
lent and had a momentum thickness of 0.0793 cm. The last measured station was at 
x = 20.95 cm or x/6 - 264. The dimensions of the cavity and of the orientations of the 
axes X,Y are not clear in reference 13. Depending on the details of the geometry, sig- 
nificant backflow and transverse pressure gradient may have been present. In view of 
these uncertainties , it is difficult to comment on the disagreement of the calculations and 
of the experiment and also on the lack of spreading of the experimental profile on the low 
velocity side. 

In any case, as the boundary layer was turbulent at the separation the shear layer 
was probably not in the self-preserving range even at the last station. (See Bradshaw, 
ref. 14.) In the calculations the maximum shear first increased above the "fully 
developed" value (overshoot). Its decrease toward the value far downstream was slow 
and monotonical. At x = 20.95 cm, its value was still higher by an amount on the order 
of 10 percent. Also the spreading rate 1/a was still substantially higher at that station 
than far downstream . 


Test Case 13 (Bradbury) 

The calculations were not overly sensitive to the initial shear. An increase of the 
initial shear by 66 percent produced less than 1-percent difference in W at x/D = 300. 

Test Case 14 (Chevray and Kovasznay) 

In a wake very close to the trailing edge, the structure of turbulence may be expected 
to be that of a boundary layer rather than that of a free flow. The question remains how 
long does it take for the structure to undergo the change from one regime to another. The 
examination of the data of reference 15 showed that the mean velocity and turbulent quan- 
tities take on their wakelike shapes already by the distance x = 20 cm and possibly even 
earlier. It seems reasonable to expect that also the structure of turbulence came close 
to the wake type by that distance. This possibility was tested by starting the calculations 
at three different points Xq = 0, 20, and 50 cm by using the experimental profiles as the 


558 



APPENDIX - Concluded 


input. The results do confirm the expectation since although the latter two calculations 
differ from the first one they agree with each other. (See figs. 12(a) and 12(b).) It 
appears that the present model has difficulties near the trailing edge where the struc- 
ture is still that of the boundary layer. The calculated and experimental profiles at 
x * 240 cm agree well when the calculations start at x = 20 cm (fig. 12(c)) or at 
x = 50 cm. From the figures it is also clear that the momentum thickness of the com- 
puted and measured profiles are not equal. A check on the momentum thickness at sev- 
eral stations revealed the following variation of 9 and gives an idea about the general 
accuracy of the experiment: 


X 

0 

20 

50 

240 

e 

0.578 

0.606 

0.585 

0.557 


Test Case 16 (Demetriades (ref. 16)) 

The results of the calculations were found to be rather sensitive to the initial shear 
in contrast to the test case 13. A 20-percent change of the initial shear produced a sub- 
stantial difference at x/40 = 2000 and this difference ^as represented by the value of 
r m /ku still persisted very far downstream. This behavior was found both in compress- 
ible and incompressible flows. This lack of tendency "to forget" of small deficit wakes 
is very interesting and suggests that convection and diffusion of turbulence in these flows 
have a large influence on the flow development. 

The presented results were obtained by using a recovery factor of 0.9. 
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Figure 8.- Small-deficit two-dimensional wake - Townsend. 






















(a) Center- line velocity. 
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(b) Maximum shear. 

Figure 12.- Two-dimensional wake of a thin flat plate 











Figure 13.- Test case 1. Free shear layers - effect of the velocity ratio. 



Figure 14.- Test case 2. Free shear layers — Mach number effect. 


570 













°l ° 


O 200 400 600 800 1000 1200 1400 1600 1600 

A 

Figure 19.- Test case 16. Two-dimensional supersonic wake. 
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DISCUSSION 


B. E. Launder : I would like to ask whether your method can be extended to cope with 
axisymmetric flows? 

T. Morel : Presumably, one could certainly try to do that, and I would expect it would 
work. However, there we have to talk, strictly speaking, about an infinite number of 
interactions. That means we would lose the nice mental picture we were basing this on. 

B. E. Launder : Doesn't it look then, that yours is a rather complicated way of doing 
something very simple ? A more conventional simulation of the shear stress equation 
would permit one to treat both axisymmetric and plane shear flows. 

T. Morel : The flows that we want to calculate are certainly not very simple and that is 
the reason we have all gathered here. We know so much more about the kinetic energy 
equation that we thought it worthwhile to pursue our work in this direction. Further, as 
I pointed out in the presentation, the interaction approach has a very important conse- 
quence. It allows us to use the kinetic energy equation to close the system without having 
to rely on the eddy viscosity to obtain the shear stress. This fact alone makes this work 
certainly worthwhile . 

G. L. Mellor : One comment here - it looks like you are taking a perfectly good energy 
equation and turning it into a shear stress equation. And yet there closely exists a per- 
fectly good shear stress equation. I think you get into trouble when you do that, as evi- 
denced by qualitative argument required to avoid the jump in sign of a, when going from 
one sign at a channel to the other. 

T. Morel : First, the exact shear stress equation is not necessarily a perfectly good 
equation. There are terms which we do not know enough about. The kinetic energy equa- 
tion is very well documented; our results certainly seem to support that. Second, the 
pressure-rate-of-strain correlation in the exact shear stress equation is usually mod- 
eled as a sum of production and dissipation. When all the terms are modeled, that equa- 
tion looks the same as our equation'. You can ask Brian Launder about that. To your last 
question, if you view the flow from the point of view of separate layers, there is no jump 
of aj within either one of them. And that is precisely the point we are making. 

A. Roshko : Yes, I would like to comment on Launder’s comment. This question of 
whether what works for two-dimensional flow will work with axisymmetric flow is not 
so clear. For example, I have a feeling that these shear layers, in particular, have the 
large structure that has a lot of two -dimensionality in it. In fact our measurements show 
that. In a fully developed axisymmetric flow, I think that is going to have a very different 
structure. I think it is an instability structure. I don't think it will be axisymmetric. I 
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think it will be skewed, and random. Therefore I am sure that the physics will be the 
same for those two flows. 

P. T. Harsha : Your figure 11 shows some pretty mystifying wiggles in your parameter 
Au/uj. Can you explain them? 

T. Morel: Well they don't mean much. You start out with some initial conditions which 
are away from self-preservation and watch what happens. K it is a flow that likes the 
self-preservation, it will tend toward it. We started a bit off, and it had to adjust. It is 
just trying to adjust. That's explainable. 

P. A. Libby : Tom, I would like to ask a question, not directly to you, but I think it does 
raise a question about some of these newer methods, and perhaps some of these other 
people will straighten me out. For example, in the method you described, if I look at 
the mathematical structure, I see a first-order differential of t with respect to y. 
That raises the possibility of satisfying one boundary condition with respect to t on 
some line of x. In a free shear problem, of course, you want to say that t is 0 at 
two points, plus or minus infinity. I don’t see how you do that. I've raised similar 
questions with other people; they say you've got to put in molecular viscosity which puts 
you back into a second-order equation. But of course, if you look at similar free-mixing 
flows, you cannot leave that molecular viscosity in, because it destroys the similarity. 
That term wants to vary as Vx and the purely turbulent problem wants to go like x 
itself. Now this, in my view, is just a manifestation of one of the problems that enter 
when you look at the newer methods of solving turbulent shear problems. I would like 
to hear what you and other people have to say about this matter. Would you like to 
comment? 

T. Morel : No, I really didn't quite follow what you said. I would like to talk to you 
afterwards. 
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